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Abstract 

In this paper we describe a simple numerical approach which allows to study the 
structure of steady-state axisymmetric relativistic jets using one-dimensional 
time-dependent simulations. It is based on the fact that for narrow jets with 
v z « c the steady-state equations of relativistic magnetohydrodynamics can be 
accurately approximated by the one-dimensional time-dependent equations after 
the substitution 2 = ct. Since only the time-dependent codes are now publicly 
available this is a valuable and efficient alternative to the development of a 
high-specialised code for the time-independent equations. The approach is also 
much cheaper and more robust compared to the relaxation method. We tested 
this technique against numerical and analytical solutions found in literature as 
well as solutions we obtained using the relaxation method and found it 
sufficiently accurate. In the process, we discovered the reason for the failure of 
the self-similar analytical model of the jet reconfinement in relatively flat 
atmospheres and elucidated the nature of radial oscillations of steady-state jets. 

Keywords: jets; relativity; magnetic fields; hydrodynamics; numerical methods 


1 Introduction 

Highly collimated flows of plasma from compact objects of stellar mass, like young 
stars, neutron stars and black holes, as well as supermassive black holes residing 
in the centers of active galaxies is a wide-spread phenomenon which has been and 
will remain the focal point of many research programs, both observational and 
theoretical. Some features of these cosmic jets, like moving knots, are best described 
using time-dependent fluid models. However, most of these jets have sufficiently 
regular global structure, which is indicative of steady production and propagation 
and promotes development of stationary models. Such models are also easier to 
analyse, and they are very helpful in our attempts to figure out the key factors of 
the jet physics. 

The simplest approach to steady-state flows is to completely ignore the variation 
of flow parameters across the jet. This allows to reduce the complicated system of 
non-linear partial differential equations (PDEs) describing the jet dynamics to a set 
of ordinary differential equations (ODEs) which can be integrated more easily [e.g. 
1, 2]. A similar reduction in the dimensionality is achieved in self-similar models, 
where unknown functions depend only on a combination of independent variables 
known as a self-similar variable. This also allows to reduce the original PDEs to a 
set of ODEs [e.g. 3, 4]. While providing important test cases and useful insights, 
this approach is not sufficiently robust - boundary and other conditions that select 
such exceptional solutions are not always present in nature. 

As it is well known to engineers working on aircraft jet engines, supersonic jets 
naturally develop quasi-periodic stationary chains of internal shocks, similar to what 
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is shown in Figure 1. These shocks emerge as a part of the adjustment of the jet 
pressure to that of the surrounding air. Interestingly, bright knots are often seen in 
cosmic jets and they are often interpreted as shocks [e.g. 5-9]. Some of these knots 
are known to be travelling and they must be part of the jet’s non-stationary dynam¬ 
ics. Others appear to be static and hence connected to the underlying quasi-steady- 
state structure of these cosmic jets. Quite often, the knots form quasi-periodic 
chains, reminiscent of those seen in aerodynamic jets. If the similarity is not ac¬ 
cidental, then these knots are also related to the process of pressure adjustment. 
In particular, we expect the powerful cosmic jets to be expanding freely soon after 
leaving their central engines and to become confined by external pressure again only 
much later [e.g. 6, 10]. The first shock driven into the jet by the external pressure is 
called the reconfinement shock. Given the growing observational evidence of station¬ 
ary knots in cosmic jets, there has been a increase of interest to the reconfinement 
process among theorists in recent years [e.g 11-17]. One of the key aims of these 
studies was to come up with approximate analytical or semi-analytical solutions for 
the structure of steady-state jets. 

Obviously, such shocked flows cannot be described by one-dimensional (ID) and 
self-similar models, which we mentioned earlier, and more complex, at least two- 
dimensional (2D), models have to be applied instead. The system of steady-state 
equations of compressible fluid dynamics, not to mention magnetohydrodynamics, 
is already very complicated and generally requires numerical treatment. One of 
the ways of finding its solutions involves integration of the original time-dependent 
equations in anticipation that if the boundary conditions are time-independent then 
the time-dependent numerical solution will naturally evolve towards a steady-state 
[e.g. 18-20]. One clear advantage of this approach is that it allows to use standard 
codes for time-dependent fluid dynamics. Such codes are now well advanced and 
widely available. However, this type of the relaxation approach is characterised by 
slow convergence and hence rather expensive. 

In order to speed up the convergence, one can use other relaxation methods, 
which are developed specifically for integrating steady-state equations [e.g. 21]. They 
often involve a relaxation variable which is called “pseudo time”. However, this 
time evolution is not realistic but designed to drive solutions towards a steady- 
state in the fastest way possible. The only disadvantage of this approach is that 
it involves development of a specialised computer code dedicated to solving only 
steady-state problems. The authors are not aware of such codes for relativistic 
hydro- and magnetohydrodynamics. 

For supersonic flows, the system of steady-state equations turns out to be hy¬ 
perbolic, with one of spatial coordinates playing the role of time [22]. (In the case 
of magnetic jets, the speed of sound is replaced with the fast magneto-sonic speed 
and we classify flows as sub-, tran-, or super-sonic based on its value compared to 
the flow speed.) In this case, one can find steady-state solutions utilising numerical 
methods which were designed specifically for hyperbolic systems, like the method of 
characteristics or “marching” schemes. These methods have been used in the past 
in applications to relativistic jets [e.g. 6, 23-26] but publicly available codes do not 
exist yet. Their development is as time-consuming as that of time-dependent codes 


Komissarov et at. 


Page 3 of 24 


whereas the range of applications is much more limited. This explains their cur¬ 
rent unavailability. Moreover, when flow becomes subsonic, even very locally, this 
approach fails. 

In this paper, we propose a new approach, which allows to find approximate 
numerical steady-state jet solutions rather cheaply and using widely available com¬ 
puter codes. To be more precise, we focus on highly relativistic narrow axisymmetric 
jets and show that in this regime the 2D steady-state equations of Special Rela¬ 
tivistic MHD (SRMHD) are well approximated by ID time-dependent equations of 
SRMHD. Like in the standard marching schemes, the spatial coordinate along the 
jet plays the role of time. This allows us to find steady-state structure of axisym¬ 
metric jets by carrying out basic ID SRMHD simulations, which can be done with 
very high resolution even on a very basic personal computer. In such simulations, 
no special effort is needed to preserve the magnetic field divergence-free and the 
computational errors associated with multi-dimensionality are eliminated. As the 
result, more extreme conditions can be tackled. Here we focus only on relativistic 
jets, because of our interest to AGN and GRB jets, but we see no reason why this 
approach cannot be applied to non-relativistic hypersonic jets as well. Our approach 
is closely related to the so-called “frozen pulse” approximation, which also utilises 
the similarity between the steady-state and time-dependent equations describing 
ultra-relativistic flows [27-29]. In this approximation, the steady-state equations 
are used to analyse the dynamics of time-dependent flows. The similarity between 
ID time-dependent models and 2D steady-state jet solutions has been noted before, 
in particular in [30]. 

In order to study the potential of this new approach we have carried out a number 
of test simulations and compared the results obtained in this way with both analyt¬ 
ical models and numerical solutions obtained with more traditional methods. The 
results are very encouraging and allow us to conclude that this method is viable 
and can be used in a wide range of astrophysical applications. 

2 Approximation 

We start by writing down the time-dependent equations of Special Relativistic Mag¬ 
netohydrodynamics (SRMHD). In this section we use units where the speed of light 
c = 1 and the factor 1/47r does not appear in the expression for the electromag¬ 
netic energy density. The components of vectors and tensors are given in normalised 
bases. The evolution equations of SRMHD include the continuity equation 


d t pT + V • (pTv) = 0 


(i) 


the Faraday equation 


d t B + VxE = 0. 


( 2 ) 


and the energy-momentum equation 


dtT tn + = o , 


(3) 
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where 



( 4 ) 


is the total stress-energy-momentum tensor, 


T h 7 = w// + pg^ 


( 5 ) 


is stress-energy-momentum tensor of matter and the components of the electromag¬ 
netic stress-energy-momentum tensor are 


T« = ( E 2 + B 2 )/2 


(6) 


T u = (ExB) 


( 7 ) 


Til = ~{E i E j + B l B j ) + E 2 + B 2 )g ij . 


(8) 


In these equations, B and E are the vectors of magnetic and electric fields respec¬ 
tively, p, p and w are the thermodynamic pressure, rest-mass density of matter and 
relativistic enthalpy of matter respectively, v is the velocity vector, T is the Lorentz 
factor and g is the metric tensor of space. These equations are to be supplemented 
with Equation of State w = w(p,p) and the Ohm’s law of ideal MHD 


( 9 ) 


E = -vxB. 


Finally, the magnetic field is divergence-free 


( 10 ) 


V-B = 0. 


In this analysis, we focus on axisymmetric jets and adopt a cylindrical coordinate 
system with the z axis coincident with the jet symmetry axis. We consider only 
narrow jets, so that 



(ii) 


We also constrain ourselves with a relatively simple magnetic configurations where 
the divergence-free condition leads to 



B r 


( 12 ) 


In axisymmetry, the steady-state Faraday equation implies E& = 0. When combined 
with Eq.9, this result yields 



v 


( 13 ) 
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Thus, the radial components of both the magnetic field and the velocity vectors are 
small compared to their axial components. 

We also assume that v^ 1. In fact, in the case of magnetically accelerated jets, 

v 4, ~ (r lc /r) 

when r r lc , the radius of light cylinder (See eq.66 in [19]). Thus, this is a good 
approximation for astrophysical jets. For a highly relativistic flow, the condition 
v z v r : v^ means 

v z ~l. (14) 

Following the standard flux freezing argument, along the jet B^/B z ~ (rj/r lc ) _1 , 
where and r s is the jet radius. (This argument does not apply to turbulent jets, which 
are non-axisymmetric and allow non-trivial conversion of components.) Hence one 
may argue that far away from the central engine 

B*^ B z . (15) 

In order to introduce the key idea of our approach we consider first the steady- 
state continuity equation: 

d z (pTv z ) + V r {pTv r ) = 0. (16) 

Using the condition (14) we may replace v z with unity. This makes Eq.16 identical 
to the ID time-dependent version of the continuity equation. In order to stress this 
point we replace 2 with t and write: 

5 t(pr) + v r (pru r ) = o. (17) 

Similarly, all 2D steady-state equations can be approximated by the corresponding 
ID time-dependent equations. 

Let us show this for the equations of magnetic field. The ID version of the diver¬ 
gence free condition reads 

d r (rB r ) =0 or rB r = const. 

Thus if B r vanishes outside of the jet, which is expected when it is in direct contact 
with ISM, then one has to put B r — 0 everywhere in the ID model. As we shell 
see, the terms involving B r are sub-dominant in all other equations and hence this 
is a reasonable simplification. Moreover, once the ID solution is found, one can 
substitute the determined B z (r, z) into the 2D divergence free condition and solve 
it for F> r (r, z). The result can then be used to verify that F> r (r, z) B z (r , z). 

The (p component of the Faraday equation can be written as 

df B 4 - rB'dj + d^B 4 ) = 0, (18) 
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where i = r,z. In steady-state, the first term vanishes, the next two terms are 
of the order B z v^/z and small compared to the last two terms, which are of the 
order B^v z /z. Removing these small terms we obtain the approximate steady-state 
equation 

d z (v z B*) + d r {v r B*)= 0. (19) 

Finally, we replace v z with unity, z with t, and obtain 

dtB* + d r {v r B*) = 0. (20) 

This is indeed the ID version of the component of eq.18. Now consider the z 
component of the Faraday equation, 

d t B z - B i diV z + -diirv'B*) = 0. (21) 

r 

The last two terms of this equation are of the order v z B z jz ~ B z jz. On the other 
hand, the second and the third terms are much smaller because of the special status 
of v z , which is approximately constant, and hence B z d z v z B z (v z /z). Removing 
these small terms, we obtain the approximate steady-state equation 

d z (v z B z ) + ld r (rv r B z ) = 0. (22) 

Now once again we replace v z with unity and z with t to obtain 

d t B z + -d r (rv r B z ) = 0, (23) 

r 

which is the ID version of the z component of eq.18. 

Finally, we analyse the energy-momentum equations. These can be written as 

d t T^ + g z T z ^ + V r T rM = 0 . (24) 

so the steady-state versions are 


d z T z * + V r T r/i = 0. (25) 

These already have the same form as the ID time-dependent equations, so we only 
need to show that 

T ZV ^ jrtu _ (26) 

Let us start with the hydrodynamic contribution. First, we notice that 

= wT 2 — p ~ wT 2 as T 1; 


Ti z d = wr2 v z - wr2 as v z ~1. 


Komissarov et at. 


Page 7 of 24 


Thus, X[f d — T* d . Then we notice that 

T£ = wTV ; 

T z d r = wT 2 v z v r ~ ioTV ; 

T h z d z = ™rVu z ~ wTV ; 

Thus, T« ~ Til 

Now we inspect the electromagnetic contributions. First, we find good estimates 
for the components of electric field. From eq.9 it follows that 

E r ~ B* (27) 

and 

E z = B r v* - 5 V <C E r . 

In fact, it is easy to show that 

E z ~ -B V . (28) 

Indeed, for magnetically accelerated jets B' l> ~ QrB z (e.g. [19]) for r 2> r lc . Hence 
v r B* ~ u r ftrB z ~ (r/r lc )B r > B r » BV . 

Using these estimates we find that 

Tt =\{E 2 + B 2 )~Bl-, 


(ExB) z ~ E r B (t> ~ B ^, 

and hence ~ T u . Moreover, 

em em / 


= -(£ 2 + s 2 ) + ^ s 2 , 

and hence T zz ~ T e tz as well. Next we show that ~ T ( U. Indeed, 

T t0 = ^ —E r B z , 

em ' 


and 


Tfjf = -{E z E <t> + B z B <t> ) ~ -B z £ r . 
Finally, we show that Tjj’ ~ TfU First, we find straight away that 


Tli = -E Z B^ and T^ =-(E z E r + B z B r ). 
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Since E z E r ~ E Z B^, we only need to show that B z B r is significantly smaller 
compared to these terms. This is indeed the case as B z B r ~ v r B z whereas using 
Eqs.27 and 28 we obtain E z E r ~ v r B^ v r B z . 

Thus, within our approximation the steady-state 2D equation of energy- 
momentum reduces to 


+ VrT r M = q ^ 


(29) 


which is the ID time-dependent energy-momentum equation. 

Given that in relativistic fluid dynamics small differences between the magnitudes 
of energy and momentum may result in huge variations of Lorentz factor and even 
lead to inconsistency, one could feel uneasy about the approximations we make. 
However, the final result is exactly the system of ID time-dependent SRMHD and 
this means that self-consistency is not compromised. For example, the flow speed 
will not exceed the speed of light because of the errors of our approximation. 

Our approach is similar to “marching” - we compute solution for a downstream jet 
cross-section using only the previously found solutions for upstream cross-sections. 
Strictly speaking, this requires the flow to be super-sonic for unmagnetised jets and 
super-fast-magnetosonic for magnetised ones [23, 31]. However, in our derivations 
we never had to utilise this condition. This suggests that it is not required when 
we wish to find only approximate solutions. For example, one may argue that the 
fact that information can propagate upstream does not necessarily imply that this 
always has a strong effect on the flow - the upstream-propagating waves could be 
rather weak. If so, we may still apply our method to jets where the supersonic 
condition is not fully satisfied, but we always need to check that the conditions 
(11-15) of our approximation hold for obtained solutions. 

3 Numerical Implementation 

The analysis of Section 2 shows that as long as they are applied to narrow jets 
with high Lorentz factor, the axisymmetric steady-state equations of SRMHD are 
very close to ID time-dependent equations of SRMHD in cylindrical geometry. This 
suggests that it may be possible to use time-dependent simulations with ID SRMHD 
codes to study the 2D structure of steady-state jet solutions. However in order to 
be able to do this, we also need to find a way of accommodating the 2D boundary 
conditions of steady-state problems in such simulations. 

For 2D supersonic flows we need to fix all flow parameters at the jet inlet and 
impose some conditions at the jet boundary, consistent with it being a stationary 
contact wave. No boundary conditions are needed for the outlet boundary - its 
flow parameters are part of the solution. In the corresponding ID problem, the 2D 
boundary conditions at the inlet boundary simply become the initial conditions 
of the ID Cauchy problem. The final ID solution corresponds to the slice of the 
2D solution at the outlet boundary. As to the contact discontinuity at the 2D jet 
boundary, the situation is not that trivial. 

Suppose that the total pressure at this boundary is a function of z, p = p h (z). 
When we replace 2 : with t this becomes p = Pb(t). Thus we need somehow to 
impose time-dependent boundary conditions. In the simulations presented below, 
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the following approach was utilised: 1) we extend the computational domain so 
that it includes the external gas, 2) we track the point separating the jet from the 
external gas and 3) we reset the external gas parameters according to the prescribed 
functions of time every computational time step. 

In order to locate the boundary separating the jet from the external gas, we 
employ a simplified version of the level-set method [32, 33]. Namely, we introduce 
the passive scalar r, which satisfies the advection equation 

d t {Tpr) + -d r (rTv r pr ) = 0 . (30) 

r 

The initial solution has a smooth distribution of this scalar 

r ~ \ (i" tanh T ~x 1 ) ’ ( 31 ) 

with the value r = 0.5 corresponding to the jet boundary (In the test simulations, 
we used A = 0.3rj.). During the simulations, the condition r < 0.5 was used to 
identify the external gas. 

After the reset, the ID jet boundary is no longer a contact but a more general 
discontinuity. In particular, the jet plasma will generally have radial velocity com¬ 
ponent. If it is positive, but in the external gas it is set to zero, then a shock wave 
will launched into the jet when this discontinuity is resolved. If it is negative, then 
this will be a rarefaction wave. On the one hand, this reflects how the information 
about changing environment is communicated to the interior of a steady-state jet. 
On the other hand, in ID simulations the strength of the emitted wave depends 
on the external density - higher density, and hence lower temperature, will result 
in stronger waves moving into the jet. This is obviously not so for 2D steady-state 
jets, which react only to the external pressure. Thus additional measures need to be 
undertaken. First, in order to negate the effect of the radial velocity jump at the jet 
boundary, the radial velocity of the external gas is reset not to zero but to its value 
at the last jet cell. Second, in order reduce the role of the external gas inertia, it 
helps to set its density to a low value, so that its sound speed becomes relativistic. 
Although we have not tried this, one could set the polytropic index of the external 
gas to T = 2, which would make the sound speed of ultra-relativistically hot gas 
equal to the speed of light. 

4 Examples 

4.1 Bowman’s jet 

To test the validity of our approach, we first use our method to reproduce the nu¬ 
merical steady-state solutions for supersonic unmagnetised jets obtained by Bow¬ 
man [25, B94] using the marching scheme described in [24]. In this study pressure- 
matched uniform jets with zero opening angle are injected into an atmosphere with 
the pressure distribution 



p(z) = Po 


(32) 
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Figure 1 Reconfinement of the Mj = 15, Tj = vTb x 10 13 K jet. The top panel is a reproduction 
of Figure 3 from B94. The bottom panel shows the solution obtained with our method. In each 
panel, the top halves show 50 pressure contours (spaced by the factor of 1.18) and the bottom 
halves show the temperature parameter r = ph/(ph — p) in 50 contours (spaced by the factor of 
1.003). The light grey lines are streamlines. 


with z s = 10, = 50. According to this equation, the external pressure initially 

decreases almost as fast as oc z~ 2 but at z > z c becomes uniform. The initial jet 
radius ro = 1 and the injection nozzle is located at z = z s . The equation of state is 
that of Synge [34] for an electron-proton plasma. The initial jet pressure pj = po. 
For the comparison we selected the model with the Mach number M 5 = 15 and 
the initial temperature T s = a/IO x 10 13 K. At such a high temperature the EOS of 
electro-proton plasma is almost the same as that of the pure proton gas. The latter 
was used in our simulations. 

Bowman’s solution is shown in the top part of figure 1. As the external pressure 
decreases rapidly, the jet quickly becomes under-expanded and enters the phase 
of almost free expansion. When it enters the outer region of constant pressure it 
becomes over-expanded and a reconfinement shock is pushed towards its axis, where 
it gets reflected. Gas passed though these two shocks becomes hot and its pressure 
rises. As a result, the jet becomes somewhat under-expanded again and begins to 
expand for the second time. Then it becomes over-expanded again and another 
shock is pushed into the jet and so on. 

In the bottom part of this figure, we show the results of our ID simulations for 
this jet using exactly the same visualisation technique as in the original paper. The 
agreement between the two solutions is quite remarkable. A very good match for 
the maximal radial extension and the oscillation-length of the jet is obtained. The 
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Figure 2 Ultra-relativistically hot jets [15, 17] in power-law atmospheres with k = 8/3, 7/3 and 2 
(from left to right). The colour-coded images show the distribution of the Lorentz factor. The 
initial Lorentz factor is Tq = 50 and opening angle Oq = 1/To- 


successive reconfinement shocks are somewhat sharper than in B94, most likely due 
to the application of a shock-capturing scheme. We checked our approach against 
other numerical models of B94 as well. In all models, the results for profile of 
jet radius and Mach number are in good agreement. Noticeable but still minor 
differences arise only for the colder models, most likely due to the different equation 
of state used in our simulations. 


4.2 Self-similar models of jet reconfinement 

The problem of reconfinement of initially free-expanding steady-state jets is quite 
important and a number of authors have tried to find simple analytic of semi- 
analytic solutions. Falle [35] and Komissarov & Falle [10] used the Kompaneets 
approximation, which assumes that the gas pressure immediately downstream of 
the reconfinement shock is equal to the external pressure at the same distance, to 
derive a simple ODE for the shock radius. Assuming particular flow profiles in the 
shocked layer, one can also determine the location of the jet boundary [e.g. 13]. 
The Kompaneets approximation is accurate only for very narrow jets. To improve 
on it, one also has to take into account the variation of the gas pressure across 
the shocked layer [11]. In our second test, we compare our results with the semi- 
analytical model by [15, thereafter KBB12], who assumed self-similarity of the flow 
in this layer. This assumption is more suitable for the case where the reconfinement 
shock never reaches the jet axis, because otherwise the distance where this occurs 
sets a characteristic length scale. 
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Figure 3 Ultra-relativistically hot jets [15, 17] in power-law atmospheres with k = 8/3, 7/3, and 
2 (from left to right). The colour-coded images show log 10 (pp _7 ). The dark blue region along 
the jet boundary is obviously a numerical artifact as its entropy is lower than that of the initial 
solution anywhere on the grid. 


KBB12 studied jets with ultra-relativistic equation of state (w = 4 p, 7 = 4/3 ), 
propagating in a power-law atmosphere, 

p=K (i) ■ (33) 

These jets emerge from a nozzle at z = zo with the Lorentz factor T 0 , opening angle 
#0 = 1/T 0 and pressure po- The initial velocity distribution correspond to a conical 
flow originating from z = 0 and hence the initial jet radius ro = zotan(#o)- They 
could only find self-consistent solutions for 8/3 < k < 4 and later argued that for 
k < 8/3 the entropy of the shocked layer must increase with the distance along the 
jet in order for the solution to be consistent with the energy conservation [17]. They 
proposed that this additional heating is caused by multiple shocks driven into the 
flow as it gradually collimates. 

We selected the KBB12 model with k = 8/3 and Tq = 50 and made simulations 
on a uniform grid with only 300 cells (each run took only several CPU minutes on a 
laptop using only one core of its processor). Our results are shown in the first panel 
of figure 2 , which should be compared with figure 7 in KBB 12 . Again we find a very 
good agreement between the models - at z = 9 we have got the jet radius r s « 0.114 
and the shock radius r s ~ 0.07, whereas in KBB12 = 0.110 and r s = 0.064. 

In order to understand the difference between the cases with k > 8/3 and k < 8/3, 
we also computed models with k = 7/3 and 2 - the evolution of the Lorentz factor 
in these models is shown in the second and third panels of figure 2 respectively. In 
these plots we see no evidence of the additional shocks proposed in [17]. Neither 
could we find them in plots of other parameters. However, figure 2 suggests that 
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Figure 4 Accuracy of the ultra-relativistic hot jets solution for the model with atmospheres with 
k — 2. The left panel shows the Lorentz factor at z — 9 for models with 150 (dotted), 300 
(dot-dashed), 600 (dashed line) and 1200 (solid line) grid points. The right panel show the total 
energy flux as a function of the distance from the nozzle for the model with 300 grid points. 


in the models with k = 7/3 and 2 the reconfinement shock is much stronger than 
in the model with k = 8/3. Moreover, the shock strength is increasing with the 
distance along the jet. As the result, the entropy of the shocked layer in the models 
with k = 7/3 and 2 is higher and its mean value across the layer is growing with the 
distance. This is confirmed in Figure 3, which shows the entropy distribution for 
these models. Since KBB12 assumed isentropy of the flow in the shocked layer, this 
could be the reason why their self-similar model fails for k < 8/3. In contrast, in 
the model with k = 8/3 the mean entropy of the layer does remain fairly constant. 
Based on these results, we conclude that the value of k = 8/3 is not special, but the 
accuracy of the constant-entropy approximation used in KBB12 greatly reduces as 
k decreases. 

The plots in Figure 3 also reveal a thin layer of decreased entropy stretching along 
the jet boundary. As in this layer the entropy is lower than anywhere in the initial 
solution, this is definitely a numerical artifact. We have checked that it becomes 
less pronounced with increased numerical resolution. Moreover, this layer forms well 
inside the jet and thus its origin is not related to the resetting procedure but is a 
property of our time-dependent code. 

We choose the model with k = 2, to illustrate the convergence and accuracy of our 
numerical solutions. The left panel of Figure 4 shows the Lorentz factor distributions 
found at z = 9 for runs with different number of grid cells in the computational 
domain, increasing from 150 to 1200 cells. As one can see, the solutions converge 
as in a first-order accurate scheme. The right panel shows the evolution of the total 
energy flux along the jet. It remains fairly constant, as expected for a conserved 
quantity. As the jet boundary jumps from one cell to another, a low level noise is 
introduced to this integral variable. 

4.3 Magnetised jets. ID versus 2D solutions 

The steady-state structure of magnetised jets is more complex, mainly due to the 
non-trivial contribution of the magnetic tension to the force balance. A number of 
authors have tackled this problem analytically using various approximations [e.g. 
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16, 36-38]. However, none of these studies deliver a model suitable for detailed 
testing of our numerical approach. Dubai & Pantano [31] studied the steady-state 
structure of relativistic jets with azimuthal magnetic field using the method of 
characteristics. This would be a good test case but the setup of their simulations 
is ambiguous. We have tried several variants of the setup but each time failed to 
reproduce the results. The mechanisms of magnetic collimation and acceleration of 
relativistic jets were studied numerically by [39, 40] and [41] using a “rigid wall” 
outer boundary. While this allows for a well-controlled experiment, Komissarov et 
al. [19] have shown that the connection between the shape of the boundary and 
the external pressure gradient is not straightforward, with significant degeneracy. 
For this reason, we concluded that in the magnetic case the best way of testing the 
performance of our ID approach would be via new 2D axisymmetric time-dependent 
simulations using the relativistic AMRVAC code [42, 43]. 

The problem we selected for this test is similar in its setup to the one described in 
Sec. 4.2 as it also involves a jet propagating through the atmosphere with the power- 
law pressure distribution (33), and the nozzle is still located at 2 = zq . However, this 
time the jet is magnetised and the rest mass density of its particles is not negligibly 
small. The jet structure at the inlet is that of a cylindrical jet in magnetostatic 
equilibrium, which satisfies the following force balance equation 

dp t b ^ drb^ 
dr r dr 

where b^ /T is the azimuthal component of the magnetic field as measured 

in the fluid frame using normalised basis and p t is the sum of the gas pressure and 
the magnetic pressure due to the axial magnetic field B z [44]. Equation (34) has 
infinitely many solutions - given a particular distribution for b^(r) one can solve 
this equation for the corresponding distribution of the pressure pt{r). We adopted 
the “core-envelope” solution of Komissarov [44] : 


6*(r) = 


/^m) ■> ^ ^ ^m 

bm(r m /r) ; r m < r < rj 
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(35) 


(36) 


where 


/3m = Wt a = 1 - 0/Pm){r m /rj) 2 , (37) 

rj is the jet radius and r m is the radius of its core (Note a typo in the expression 
for a in [44].). As one can see, the core is pinched and in the envelope the magnetic 
field is force-free. This may be combined with any distribution of density and axial 
velocity. We imposed p = po and 


r(r) = To (1 - (r/rjY) + (r/ rj ) 


(38) 
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Figure 5 Initial radial structure of the magnetised jets in the test simulations described in Sec. 4.3. 


with v = 8 ; this gives an almost “top-hat” velocity profile. The velocity vector is 
set to be aligned with the jet axis, so v r m = 0. This solution is illustrated in 
Figure 5. 

We considered two models, A and B. In the models A, the magnetic field is purely 
azimuthal and the other parameters are 77 = 1, r m = 0.37, b m = 1, po = 1 5 
zo = 1, Pm = 0.34, To = 10. The local magnetisation parameter a = b 2 /w does not 
exceed cr max = 0.7 in this model and thus the jet is only moderately magnetised. 
The jet core is relativistically hot, with the gas pressure reaching p m ax = p at the 
axis, which opens the possibility of efficient hydrodynamic acceleration once the jet 
is allowed to expand. In the simulations we used the adiabatic equation of state 
w = p -\- (yy/'y — 1 )p with 7 = 4/3. 

In model B, this configuration is modified to include nonvanishing longitudinal 
magnetic field B z . In particular, we considered the case where the gas pressure 
p = apo everywhere within the jet and 


B z 


Po 


£(1 ~{r/rm) 2 ) 
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rri 

m •> 


(39) 


which keeps p t unchanged. In this model, the magnetic field is force-free not only in 
the envelope but also in the core. The other parameters of this model that differ from 
those of model A are po = 0.05 and /3 m = 0.14. The corresponding magnetisation 
is much higher, with cr max = 17. 

Model B turned out too stiff for our 2D code, but presented no problems in ID 
simulations. For this reason we compare here the ID and 2 D results for model A 
only. In these simulations we used the atmosphere with k = 1. The computational 
domain is 20 77 in the radial direction and 800 77 in the axial direction. 
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First, let us describe the overall jet structure found in these simulations. Initially, 
as the jet enters the region of rapidly declining external pressure, it expands rapidly 
and a rarefaction wave moves towards its axis. Eventually, the jet becomes over- 
expanded, its expansion slows down, and a reconfinement shock sets in. It reaches 
the axis at z « 400, gets reflected and then returns to the jet boundary at z ~ 700 
(see Figure 6). 

To quantify the convergence of the ID simulations we carried out simulations with 
three different resolution and used this data to determine the grid-convergence index 



(40) 


where /i, are solutions with doubled and quadrupled resolution compared to /o 
and | f a — f\) |i is the difference between two solutions in the L\ norm. We found 
that rj « 1, as this is expected for a TVD scheme in the presence of discontinuities. 
At 6400 grid cells in the radial direction, the density contours become visually 
unchanged on the scale of figure 6. The ID solution with this resolution was used 
for comparison with the results of our 2D simulations. In what follows we refer to 
it as the “converged” ID solution. 

The initial solution in our 2D simulations was constructed via interpolation of the 
converged ID solution onto the 2D cylindrical grid. Since we did not include gravity 
to balance the pressure gradient in the external atmosphere, in order to preserve 
the atmosphere in its initial state the atmospheric parameters were reset to their 
initial values every time step, just like this was done in the ID case. In order to 
test the convergence of 2D solutions, we made three runs with doubled resolution, 
N r = 400, 800, and 1600 cells in the radial direction. The number of cell in the 
axial direction was always twice the number of cells in the radial one. 

Typically, the 2D solutions exhibited some evolution at first but then quickly 
settled into a stationary state. For example in the case of N r = 400, the timestep-to- 
timestep relative variation of the conserved flow variables dropped below 6 x 10 -6 
at t = 1000 and remained approximately constant thereafter. Furthermore, the 
relative L\ error of density between times t = 1000 and t = 3000 was 2.8 x 10 -4 , 
indicating that a stationary state had been reached. The 2D solutions converge with 
the grid-convergence index p > 1.25 over the entire simulated time. 

The difference between the converged ID solution and the relaxed 2D solutions 
with with N r = 800 (dotted lines) and N r = 1600 (dashed lines) is illustrated 
in Figure 6 which shows the mass density distribution. One can see that the 2D 
solutions are very close to the ID solution and that the difference decreases with the 
resolution of 2D runs. To quantify the difference between the relaxed 2D solutions 
and the converged ID approximate solution we introduce the parameter 


Sp = \p2D ~ Pid\i/(pid) • 


(41) 


For the 2D solution with N r = 400 cells in the radial direction we obtain Sp ~ 6%, 
for N r = 800 Sp ~ 4.3% and for N r = 1600 the relative error decreased to Sp ~ 3.2%. 
This shows that the approximation error of our ID approach is at the level of no 
more than 3%. 
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Figure 6 Converged ID solution for a stationary magnetised jet (solid lines) and two 
corresponding 2D solutions found via the relaxation method, one with the 1600 x 3200-resolution 
(dashed lines) and one with the 800 x 1600-resolution (dotted lines). The lines are 10 rest-frame 
density contours consecutively spaced by the factor of one half from the starting value of 
Pmax = 1. The solution involves a reconfinement shock which reaches the jet axis at z ~ 400. 


4.4 Magnetised jets in power-law atmospheres 

Komissarov et al. [19] derived an approximate equation for the radius of highly 
magnetised jets, in the limit where it strongly exceeds that of the light cylinder. 
Using this equation they concluded that in the case of power-law atmosphere with 
0 < k < 2 the jet radius increases as 


,k/4 


(42) 


rj oc z 


Lyubarsky [37] developed the theory of Poynting-dominated jets further and using 
more accurate analysis found that the expansion is modulated by oscillations with 
the wavelength growing as 


A oc z K / 2 


(43) 


These oscillations can be understood as a standing magneto-sonic wave bouncing 
across the jet. Indeed, denote the wave speed as a m . Then the jet crossing time 
is r c = Vj/cim in the co-moving jet frame and t c = Tr c in the rest frame of the 
atmosphere. As the wave is advected along the jet almost at the speed of light the 
wavelength of the associated structure is 



( 44 ) 
























Komissarov et at. 


Page 18 of 24 



Figure 7 Structure of steady-state magnetised jets obtained via time-dependent ID simulations. 
The plots show the co-moving density distribution for model A with k = 0.5 and k, = 1. The 
distance along the vertical axis is defined as z = ct/rj, where rj is the initial jet radius. The white 
contour shows the jet boundary, located using the passive scalar. 


Since for the jets considered in [19, 37] a m ~ c and T oc 77 we obtain A oc r 2 and 
using Eq.42 recover Eq.43. The results (42,43) are well suited for testing of our 
approach. To this aim, we carried out additional ID simulations with models A and 
B described in Sec. 4.3. 

Since in model A the jet is not Poynting-dominated, it allows us to explore the 
regime not covered in [37]. To see how sensitive these results may be to the assump¬ 
tions made in [19, 37] let us consider unmagnetised relativistic jets. From the mass 
conservation law we obtain 77 oc (Tp)~ 2 . For relativistically cold jets with p <C pc 2 
we have T ~ const and thus 

77 oc z K t 21 , (45) 

whereas for the hot jets T oc 77 and thus 

rj oc W 4 , (46) 

where we put 7 = 4/3. The last result is the same as for the Poynting-dominated 
jets. Even for the cold jets the difference is rather minor, e.g. for 7 = 5/3 the index 
in Eq.45 differs from /c/4 only by /c/20 and for 7 = 4/3 by /c/ 8 . 

In order find A we note that for cold jets oc (p/p) oc z _/ d 1 _ 1 / 7 ) and hence 
Eq.44 yields Eq.43 independently of the value of 7 . For hot jets, a m ~ const and 
Eq.44 still leads to Eq.43 if we use 7 = 4/3. Thus, the law (43) for the wavelength 
of oscillations is very robust. 
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Figure 8 Structure of the steady-state magnetised jet in model A with k, = 1, obtained via 
time-dependent ID simulations. From left to right, the plots show the total pressure, Lorentz 
factor and the magnetisation parameter cr. Jet oscillations cause compression in the squeezed 
regions as well as re-acceleration of the bulk flow as the flow expands. The majority of 
acceleration occurs in the thermally dominated core. A reconfinement shock is clearly visible in 
the total pressure and magnetisation plots. 


Figure 7 illustrates the overall jet structure in model A and its response to changes 
in the parameter k of the external atmosphere. One can see that this weakly mag¬ 
netised jet also shows a combination of secular expansion and oscillations. These 
oscillations appear to be a generic feature of the adjustment process of supersonic 
jets to variations of external pressure, which occurs by means of magneto-sonic 
waves travelling across the jet. In the very beginning, the decrease of external pres¬ 
sure makes the jet under-expanded and a rarefaction wave is launched from the 
jet boundary towards the jet axis. Behind this wave the radial velocity is positive 
and the flow expands. The rarefaction reduces the jet pressure and at some point 
it becomes over-expanded. Now a compression wave is driven inside the jet. Behind 
this wave the jet expansion slows down and eventually turns into a contraction. The 
contraction increases the jet pressure and at some point it becomes under-expanded 
again and then the whole cycle repeats. 

The deviation from the force-balance corresponding to the secular jet expansion 
is due to the finite propagation speed of the waves - as they move across the jet they 
are also advected downstream by the supersonic flow. As the result, the jet interior 
reacts to the changes in the external pressure with a delay. It keeps expanding when 
the internal pressure is already too low and keeps contracting when it is already 
too high. As k increases, the wavelength of the oscillation increases as well. This 
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Figure 9 As figure 8 for model B with k = 1. As the Poynting-flux vanishes on the axis (and the 
thermal energy is negligible), we obtain a hollow jet with fastest region away from the axis. Due 
to the increased fast-magneto-sonic speed (thus lower Mach-number) compared to the case of 
figure 8, no reconfinement-shock forms and the jet-oscillation frequency is increased. 


is expected as the more rapid overall expansion of the jet in an atmosphere with 
larger k means that it takes longer for a magneto-sonic wave to traverse the jet, not 
only as the result of the larger jet radius but also as the result of its higher Mach 
number (and hence smaller Mach angle). 

Overall, this is very similar to the well-known evolution of under-expanded su¬ 
personic jets studied in laboratories. Normally, their compressive transverse waves 
steepen into shocks. In our model A with k = 1 we also detect shocks, but they 
become progressively weaker, suggesting that they may disappear further out along 
the jet. For k = 0.5, shocks do not form at all. The exact reason for this in not yet 
clear. 

Figure 8 shows the evolution of other flow parameters in model A with k = 1. 
Both the secular and oscillatory behaviour of the jet radius are mirrored in the 
variation of the Lorentz factor. The secular expansion leads to secular increase of 
the Lorentz factor as both the thermal and the magnetic energy are converted into 
the kinetic energy of the flow. The thermal acceleration is most pronounced in the 
jet core, which is relativistically hot at the inlet. The oscillations of the jet radius 
lead to additional increase of the Lorentz factor during the expansion phase and 
its decrease upon contraction. The left panel of Figure 10 shows the dynamics of 
energy fluxes for this jet. These are found via integration over the jet cross-section of 
b 2 V 2 v z — b°b z for the magnetic energy, pT 2 v z for the kinetic energy and (w — p)Y 2 v z 
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for the thermal energy. The results are normalised to the rest-mass flux, obtained 
via integration of pTv z . As the result of this normalisation, the kinetic energy flux 
has the meaning of mean actual Lorentz factor of the jet, whereas for the thermal 
energy and magnetic energies these are the gains in the Lorentz factor, which can be 
achieved upon full conversion of these energies into the kinetic one. The main feature 
of the plot is a conversion of the thermal energy into the kinetic one (the magnetic 
energy is highly sub-dominant from the start). This conversion is largely completed 
during the initial phase of monotonic expansion, which lasts up to z = 200. In the 
second phase, the thermal energy flux is comparable to the magnetic, and they are 
being converted to the kinetic energy at more of less the same and rather slow rate. 
Strong oscillations are superimposed upon this secular evolution, with the kinetic 
(thermal) energy reaching local maxima (minima) at the locations of jet bulging. 

Figure 9 shows the same parameters for the highly magnetised jet of model B 
with k = 1. In this model, the reconfinement shock is no longer present. This may 
be related to the fact that in this model the fast magneto-sonic speed is higher 
and the corresponding jet Mach number is lower, at the inlet M ~ 3 compared to 
M ~ 10 in model A. The lower Mach number is also responsible for the observed 
lower wavelength of the jet oscillations as it takes less time for the waves to traverse 
the jet. In this model, the jet is magnetically-dominated and the main feature of its 
energy balance is a gradual conversion of the magnetic energy into the kinetic one 
(see the right panel of Figure 10). 

The theoretical predictions for the secular evolution of the jet radius and the 
wavelength of its oscillations are put to a quantitative test in figure 11 , which shows 
the jet radius rescaled according to its expected secular evolution against . 

In such plots, the mean jet radius and the wavelength of oscillations should remain 
constant. For the highly magnetised jet of model B the scaling factor is W 4 and for 
the low magnetised jet of model A it is z 3 ^/ 8 , as appropriate for a cold hydrodynamic 
jet with 7 = 4/3. In general, we obtain a very good agreement with the theoretical 
scalings for the mean jet radius, both in the low- and high-magnetisation limit. A 
small departure from the z^-scaling is observed for case B with k = 1 - it expands 
slightly faster. This could be because the jet magnetisation is not sufficiently high 
and decreases more rapidly with distance than in the atmosphere with k = 0.5. 


























Komissarov et al. 


Page 22 of 24 


50 


100 


150 

7 1—k/2 


200 



250 


300 


Figure 11 Compensated jet-expansion laws for models A (top) and B (bottom). In both models 
the expected average expansion is captured quite well. To show that the oscillation wavelength 
scales as A oc z K ! 2 , the x-axis has been rescaled accordingly. In order to visually separate the 
curves corresponding to different values of k, they have been shifted up by a factor of k. 


The evolution of the wavelength scaling is also in a very good agreement with the 
theory - the residual error is between 0.7% and 3.4%. 

5 Conclusions 

In this paper we presented a novel numerical approach, which can be used to de¬ 
termine the structure of steady-state relativistic jets. It is based on the similar¬ 
ity between the two-dimensional steady-state equations and the one-dimensional 
time-dependent equations of SRMHD with the cylindrical symmetry in problems 
involving narrow highly-relativistic (v z ~ c) flows. Such similarity has already been 
utilised in the so-called “frozen pulse” approximation where dynamics of time- 
dependent relativistic flows is analysed using the steady-state equations [27-29]. 
Here we do the opposite and construct approximate steady-state solutions via nu¬ 
merical integration of the time-dependent equations. The main advantage of this 
approach is utilitarian. First, it allows us to use computer codes for relativistic 
MHD (or hydrodynamics in the case of unmagnetised flows), which are now widely 
available, in place of highly-specialised codes for integrating steady-state equations, 
which are not openly available at the moment. Moreover, the reduced dimensional¬ 
ity means that the computational facilities can be very modest - a basic laptop will 
suffice. In contrast, the relaxation method based on integration of two-dimensional 
time-dependent equations can be computationally quite expensive. 

We compared numerical solutions obtained with this approach with analytical 
models and numerical solutions obtained with other techniques. The considered 
problems involved a variety of flows both magnetised and unmagnetised, with dif¬ 
ferent equations of state and external conditions. The results show that the method 
is sufficiently accurate and robust. 

Although we focused only on relativistic flows, we see no reason why this approach 
cannot be applied to non-relativistic hypersonic flows. For such flows, the axial 
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velocity of bulk motion plays the role of the speed of light in the substitution z = ct 
used in our derivations. 

As a byproduct of our test simulations, we obtained two results of astrophysical 
interest. We demonstrated that the failure of the self-similar model of the jet re¬ 
confinement in power-law atmospheres with the index k < 8/3 [15] is rooted in the 
assumption of isentropy of the shocked layer, which is made in this model. In reality, 
the reconfinement shock becomes stronger with the distance along the jet, resulting 
in a strong spatial variation of the entropy. We also found that the radial oscillations 
of steady-state jets, discovered in the analytical models of Poynting-dominated jets 
[37] is a generic part of the jet adjustment to the space-variable external pressure 
and not specific to the high-magnetisation regime only. The oscillations are standing 
waves induced by the variation. 

The steady-state solutions are useful for elucidating some key factors in flow 
dynamics and may closely describe some of the observed phenomena in astrophysical 
jets. However, they are often subject to various instabilities which may dramatically 
modify the flow properties. Most instability studies, both analytical and numerical, 
deal with very simple problems where the steady-state solution is readily available. 
In more realistic setup, the issue of finding the steady-state solution, which can 
then be subjected to perturbations, becomes more involved and this is where our 
method can be applied in the instability studies. 
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